Inferring predator–prey interactions from camera traps: A Bayesian co‐abundance modeling approach

Abstract Predator–prey dynamics are a fundamental part of ecology, but directly studying interactions has proven difficult. The proliferation of camera trapping has enabled the collection of large datasets on wildlife, but researchers face hurdles inferring interactions from observational data. Recent advances in hierarchical co‐abundance models infer species interactions while accounting for two species' detection probabilities, shared responses to environmental covariates, and propagate uncertainty throughout the entire modeling process. However, current approaches remain unsuitable for interacting species whose natural densities differ by an order of magnitude and have contrasting detection probabilities, such as predator–prey interactions, which introduce zero inflation and overdispersion in count histories. Here, we developed a Bayesian hierarchical N‐mixture co‐abundance model that is suitable for inferring predator–prey interactions. We accounted for excessive zeros in count histories using an informed zero‐inflated Poisson distribution in the abundance formula and accounted for overdispersion in count histories by including a random effect per sampling unit and sampling occasion in the detection probability formula. We demonstrate that models with these modifications outperform alternative approaches, improve model goodness‐of‐fit, and overcome parameter convergence failures. We highlight its utility using 20 camera trapping datasets from 10 tropical forest landscapes in Southeast Asia and estimate four predator–prey relationships between tigers, clouded leopards, and muntjac and sambar deer. Tigers had a negative effect on muntjac abundance, providing support for top‐down regulation, while clouded leopards had a positive effect on muntjac and sambar deer, likely driven by shared responses to unmodelled covariates like hunting. This Bayesian co‐abundance modeling approach to quantify predator–prey relationships is widely applicable across species, ecosystems, and sampling approaches and may be useful in forecasting cascading impacts following widespread predator declines. Taken together, this approach facilitates a nuanced and mechanistic understanding of food‐web ecology.


| INTRODUC TI ON
Understanding predator-prey interactions is a foundational theme in ecology (Lotka, 1920;Volterra, 1927). The importance of predator-prey interactions in conservation biology has gained widespread attention following the global decline of apex predators and subsequent trophic cascades (Estes et al., 2011;Ripple et al., 2014). Predators can play keystone roles in structuring ecosystems by regulating prey populations and the spatiotemporal distribution of prey via consumptive and behavioral effects (Estes & Palmisano, 1974;Gaynor et al., 2019). Where predators exert strong top-down control and suppress prey, predatory interactions can produce a negative co-abundance relationship (Ripple et al., 2014).
For example, if prey are primarily bottom-up regulated by limited food availability, and predators are consequentially bottom-up limited by prey availability, a positive co-abundance relationship may arise, and this has been observed previously (Figure 1a; Carbone & Gittleman, 2002;Karanth et al., 2004). A positive predator-prey coabundance relationship may also suggest alternative forces may be more important than predation in regulating species' abundances (e.g., bottom-up control, hunting, or other severe disturbances ;Ford & Goheen, 2015). Despite the importance of predatory interactions in food-web ecology and conservation, measuring predator-prey relationships in natural settings imposes several key obstacles due to direct observations, diet analyses, or manipulative experiments being logistically difficult (Smith et al., 2020), especially for cryptic species in tropical forests (Brodie & Giordano, 2013).
Ecologists have frequently used observational methods to study predation by comparing landscapes that vary in predator and prey densities (e.g., Atkins et al., 2019;Brashares et al., 2010;Ripple & Beschta, 2006;Terborgh et al., 2001). However, modern sampling methods that rely on passively observed occurrence data (e.g., from camera traps, acoustic monitors, or eDNA) face numerous analytical barriers to inferring interactions (Blanchet et al., 2020;Ford & Goheen, 2015). Observational studies rarely account for processes influencing both predator and prey, such as environmental and anthropogenic covariates, the complete absence of either species from specific landscapes (i.e., true zeros), and imperfect detection of each species (i.e., false zeros), all of which can bias estimates of interactions (Blanchet et al., 2020;Blasco-Moreno et al., 2019). For example, generalized joint attribute models can accommodate zeroinflated and over-dispersed data to infer species interactions (Clark et al., 2017), but such approaches fail to account for imperfect detection. Recent attempts to use detection-corrected abundance or co-occurrence models have failed to propagate uncertainly throughout the modeling process, which may spuriously increase their statistical power (e.g., Penjor et al., 2022). This leaves a gap in our ability to routinely study predator-prey interactions using observational data from standardized biodiversity monitoring programs that now often use camera trapping (Jansen et al., 2014). Here we leverage recent advances in collecting large standardized multispecies datasets via camera trapping and hierarchical N-mixture co-abundance models to enhance our understanding of complex wildlife interactions in natural settings.
Many ecological sampling methods produce multispecies detection histories (i.e., binary detections or nondetections per sampling location and sampling occasion) or count histories (i.e., number of individuals counted per sampling location and sampling occasion), including point counts, camera trapping, and acoustic sampling.
Detection and count histories can capture the spatiotemporal patterns of species observed across a landscape and have been used

T A X O N O M Y C L A S S I F I C A T I O N
Community ecology, Trophic interactions F I G U R E 1 Hypothetical predator-prey co-abundance relationships. (a) A positive predator-prey relationship suggests top-down regulation via predation is not the dominant factor shaping abundances. This may arise where there is strong bottomup regulation of both prey and predators or severe disturbances affecting both species. (b) The lack of a predator-prey relationship could arise due to a lack of interactions, such as if predator dietary preferences exclude a specific prey species or each species utilizes different habitats. (c) A negative predator-prey relationship could arise from strong top-down regulation via predation that regulates prey abundance. Similarly, a negative predator-prey relationship suggests predator extirpation may allow prey to increase, termed "trophic release".
to study habitat associations, species distributions, and population dynamics (Kéry & Royle, 2016) but have rarely been used to measure predator-prey relationships (Kéry & Royle, 2021). Detection histories have commonly been used to examine species interactions through several varieties of co-occurrence models that may quantify nonindependent occurrence (i.e., symmetrical interactions) within the entire community (Tobler et al., 2019) or quantify a directional effect where the occurrence of one species is conditional upon the presence of another (i.e., asymmetrical interactions; Richmond et al., 2010), and both approaches account for imperfect detection.
By contrast, count histories have received considerably less attention for inferring species interactions despite carrying more information than detection histories (Roth et al., 2016). Recent work by Blanchet et al. (2020) outlined numerous concerns about using detection histories to infer species interactions and argued that occurrence data does not carry enough information to infer interactions, species occurrence depends on their environment (i.e., habitat filtering), and that strong interactions may lead to exclusions before interactions can be detected.
Co-abundance models that use count data instead of occurrence data to infer species interactions have solved many of the issues raised by Blanchet et al. (2020;Table 1;Brodie et al., 2018;Roth et al., 2016). Co-abundance models are an extension of widely used N-mixture models used to estimate abundance for species that cannot be individually identified (Royle, 2004). N-mixture models can generate latent population size estimates that are TA B L E 1 Problems and solutions for inferring species interactions from observational data, such as camera trapping capture histories.

Problems inferring predator-prey interactions Solutions
Occurrence data (i.e., presence and absence) lacks enough information to infer biotic interactions (Blanchet et al., 2020;Roth et al., 2016) Camera trapping can obtain count data (i.e., number of individuals per independent capture) and N-mixture models utilize count data in a count history matrix to estimate abundance while accounting for imperfect detection (Royle, 2004) Need to propagate uncertainty throughout the entire modeling process (Brodie & Giordano, 2013) Use integrated co-abundance model where all parameters are estimated in a single set of simulations using Bayesian Markov chain Mote Carlo (MCMC) methods (Brodie et al., 2018;Roth et al., 2016) Species occurrence and abundance depend on the environment (i.e., habitat filtering) and the joint habitat preferences of two species could falsely create the illusion of interactions (Blanchet et al., 2020;Dormann et al., 2018) Integrate environmental covariates affecting both species into the co-abundance modeling framework (Brodie et al., 2018;Roth et al., 2016) Sampling scale influences measures of co-occurrence. Spatial scale must be fine scale enough to detect interactions between individuals, while also encompassing broad variation in species' distributions (Blanchet et al., 2020) Spatially resample camera trap locations to reflect the study species' home ranges that ensure detections are spatially independent and comparable across camera trapping sessions (Rayan & Linkie, 2020). Also, use a large multi-landscape dataset that represents diverse samples of both species' distributions Appropriate statistical inference requires a very large sample size (Blanchet et al., 2020) Technological and cost improvements have made it possible to conduct large, standardized, and repeated camera trapping sessions that can collectively generate sufficient sample sizes, even for rare and cryptic species. For example, in our case study below, we used data from 1210 camera traps from 20 sessions conducted across 10 landscapes that generated 5980 independent captures of our four study species Strong interactions may lead to complete species exclusions and adding zeros to count histories usually leads to deviations from the normal Poisson distribution (Blanchet et al., 2020) Incorporate species absences in count histories and classify zeros as true or false zeros (Blasco-Moreno et al., 2019). Account for true zeros (e.g., species extirpated from a landscape) using an informed zero-inflated Poisson distribution, and account for false zeros (e.g., present but not detected) by correcting for imperfect detection Positive correlations between predator and prey abundances are difficult to interpret in terms of interspecific interactions (Brodie et al., 2018) Interpret negative predator-prey relationships as evidence that predation is regulating the focal prey population(s) (Ripple et al., 2014). Interpret positive predator-prey relationships as evidence that predation is not the primary force regulating the focal prey population(s). Positive predator-prey relationships suggest alternative forces (e.g., bottom-up control, hunting) may be more important in shaping the focal prey species abundance (Ford & Goheen, 2015) N-mixture models are sensitive to model assumption violations that can inflate absolute density estimates leading to inaccurate inferences for determining population size (Link et al., 2018) Interpret results from co-abundance models as the directional change in species abundance relative to covariates and not as absolute population sizes (Gilbert et al., 2021) Note: This manuscript describes an approach to implement all solutions.
roughly equivalent to capture-recapture analyses under ideal sampling conditions (Ficetola et al., 2018), though these models can be sensitive to assumption violations that may inflate true population density estimates (Link et al., 2018;Nakashima, 2020). The key advantage of using N-mixture models is that they accurately quantify the spatial variation in abundance as a function of covariates, thus producing a relative, rather than absolute, a measure of abundance (Gilbert et al., 2021). Currently, co-abundance models that rely on N-mixture models to infer species interactions have been limited to competing species that both occur at similar densities so both species' count histories can be assumed to follow a normal Poisson distribution (Brodie et al., 2018;Cosentino et al., 2019;Easter et al., 2020;Roth et al., 2016). Existing co-abundance models exhibit poor performance and parameter convergence when species count histories have vastly different distributions, limiting their applicability for predators and prey whose natural densities vary by an order of magnitude and exhibit considerably different detection probabilities (Carbone & Gittleman, 2002;Sollmann et al., 2013). We developed a Bayesian hierarchical N-mixture co-abundance model to study predator-prey interactions while conforming to the criteria proposed by Blanchet et al. (2020) for robust inferences (Table 1). Our key improvements over existing co-abundance models for competition (Brodie et al., 2018;Cosentino et al., 2019;Easter et al., 2020) are introducing an informed zero-inflated Poisson distribution in the abundance formula and accounting for overdispersion in detections using a random effect per sampling unit and sampling occasion in the detection formula. We illustrate that the combination of these two methodological advancements is necessary to infer ecologically meaningful predator-prey interactions while ensuring parameter convergence and goodness-of-fit across several different species pairs. Our co-abundance models facilitate the study of predator-prey interactions across trophic levels by quantifying predator-prey co-abundance relationships ( Figure 1), while accounting for each species' detection probability, relationships with environmental covariates, and propagates uncertainty. As a real example, we examined four potential predatory interactions using a multi-session multi-landscape camera trapping dataset from Southeast Asian tropical forests. We quantified predator-prey co-abundance relationships to test whether tigers (Panthera tigris) or clouded leopards (Neofelis nebulosa and N. diardi) suppress sambar deer (Rusa unicolor) or muntjac deer (Muntiacus muntjac).

| A two-species N-mixture model
Our approach was to test for the effect of a dominant predator's local abundance on a subordinate prey's local abundance using count histories. We adapted the N-mixture co-abundance modeling framework from Brodie et al. (2018) that included a term for the latent abundance of a dominant species affecting the abundance of a subordinate species. This N-mixture model estimates local abundance for species i (either dominant dom or subordinate sub ) at sampling unit j, denoted as N i,j , through repeated counts of the population over a time frame during which the population is closed to change (Royle, 2004). We assume that: and model the expected count of species i at sampling unit j, λ i,j , relative to covariates using a log-link function (Royle, 2004). However, we Therefore, our iZIP N-mixture model assumes that: Fixing local abundance to zero rather than estimating nonzero abundance where the species was extirpated minimizes the chances of making a type I error (Martin et al., 2005).
The effect of the covariates and the dominant on the subordinate species was modeled as: where α sub is a vector of environmental covariate effects for the subordinate and δ represents the effect of the latent abundance per sampling unit of the dominant species (N dom,j ) on the subordinate.
An estimated value of δ < 0 would infer a negative co-abundance relationship between the dominant and subordinate species in support of a predatory interaction (e.g., top-down regulation; Figure 1c).
An estimated value of δ > 0 would infer a positive predator-prey co-abundance relationship that is likely driven by responses to unmodelled covariates (e.g., food availability to prey that underlies bottom-up regulation, Figure 1a). Estimates of δ with 95% Bayesian credible intervals (CI) overlapping zero do not indicate biologically meaningful interactions ( Figure 1b). Assuming the covariates included in this model are appropriate, this approach allows us to tease apart the directional impact of a dominant on subordinate species (Brodie et al., 2018). The abundance model for both dominant and subordinate species included the iZIP parameter and the same environmental covariates, while the subordinates included the additional parameter estimating the effect of the latent abundance of the dom- We assume that we imperfectly observed the latent abundance of both species during sampling, thus giving rise to false zeros in our count histories (Blasco-Moreno et al., 2019;Royle, 2004).
Abundance cannot be directly observed at sampling units, so sampling biases like imperfect detection are accommodated in estimates of N i,j by assuming that the detections of species i at sampling unit j during sampling occasion k, denoted as n i,j,k , follow a binomial distribution with species-level detection probability p i,j,k : We modeled the detection probability of species i at sampling unit j during sampling occasion k, p i,j,k , relative to covariates using a logitlink function (Royle, 2004), and we used the same sampling-related covariates in the detection model for both species. The original Nmixture model with a binomial detection process assumes that individuals are not double-counted between sampling occasions and sampling units (Royle, 2004), and this assumption may be violated if the same individual is detected multiple times within the same sampling occasion. Double counting individuals can lead to inflated abundance estimates, so there have been suggestions to address this using a Poisson detection process, which we also tested (Link et al., 2018;Nakashima, 2020). To account for overdispersion in species-specific detection probability not captured by samplingrelated covariates, we included a random effect per sampling unit and sampling occasion (an overdispersion random effect, hereafter "ODRE") per species, ε i,j,k (Kéry & Royle, 2016;Roth et al., 2016). For the ODRE, we assume that: where τ i was the standard deviation of the ODRE. We also included a stabilizing parameter to ensure the logit-link transformation would not become zero or negative (Kéry & Royle, 2016). Bringing the samplingrelated covariates and ODRE together, we define the detection probability of species i at sampling unit j on sampling occasion k as: We estimated model parameters using a Bayesian approach with MCMC methods with the program R (R Development Core Team, 2021) in the package "jagsUI" (Kellner, 2019). We used this Bayesian approach to propagate uncertainty throughout the modeling process by first estimating the latent abundance of the dominant species per sampling unit, N dom,j , which is then used as a covariate to estimate the latent abundance of the subordinate per sampling unit, N sub,j (Brodie et al., 2018;Roth et al., 2016). As our primary goal for this N-mixture co-abundance model was the examine the directional impact of a dominant species' local abundance upon a subordinate species' local abundance, and not to estimate latent population sizes, we refrain from interpreting our results in terms of absolute density but rather as the spatial variation in abundance relative to covariates (i.e., δ * N domj,k ; Gilbert et al., 2021). Apart from our iZIP parameter (Z ij ), we used uninformative prior values and provided similar initial values close to zero for all parameters. For each species pair, we ran three chains of 1,000,000 iterations because Kéry and Royle (2016) stressed the importance of running long chains when incorporating the ODRE.
We discarded the first 200,000 iterations as burn-in and thinned the chains by 80, which left 30,000 values to quantify the posterior distribution of each parameter. We assumed parameters converged if their Rhat scores were between 1 and 1.2 (Gelman et al., 2013). We calculated the 95% CI from the posterior distribution and considered our parameters to have a clear effect (or "significant effect" in frequentist terminology) if 0 was not included in the 95% CI (Roth et al., 2016).
To infer confidence about parameter directionality, we calculated the probability the posterior distribution of the species interaction parameter was either negative or positive using the R package "bayestestR" (Makowski et al., 2019). Finally, our model has the same assumptions as a standard N-mixture model, including independence among sampling units, population closure over all replicated sampling occasions, independent and equal detection probabilities among individuals within a species, and that abundance per sampling unit was a random variable with E(N j ) = λ (Royle, 2004).

| Camera trapping methods
We assessed predator-prey co-abundance relationships using a large multi-session multi-landscape camera trap dataset from Southeast Asian tropical forests. We conducted 20 camera trapping sessions in ten lowland primary tropical forest landscapes in Sumatra (3), Borneo (2), Singapore (1), Peninsular Malaysia (2), and Thailand (2) (Phillipps & Phillipps, 2016). We only included count data from landscapes where the predator species is native, which allowed us to incorporate landscapes where the species is extirpated (e.g., tigers in Singapore) and exclude landscapes where the species is not native (e.g., tigers in Borneo, Amir et al., 2022). Sambar deer (Rusa unicolor) and muntjac deer (Muntiacus muntjak) deer are common large ungulates widely distributed through our study region. Petersen et al., 2020). Therefore, we expected strong negative co-abundance relationships between tigers and sambar deer and between clouded leopards and muntjac deer in support of predatory interactions.

| Covariates
While there are many factors that affect species abundance, our analysis focused on a limited set of powerful composite variables.
To account for environmental and anthropogenic disturbances that may affect species abundance, we included the Forest Landscape Integrity Index (FLII) and Human Footprint Index (HFP) as fixed effect covariates on the abundance formula for both species. The FLII is a globally continuous measure of the world's forests that integrates both observed deforestation and inferred edge effects and the loss of connectivity (Grantham et al., 2020). The HFP is a globally continuous measure that represents landscape-level anthropogenic disturbances from human population densities and infrastructure, and can be used as a crude metric to infer potential hunting pressure (Venter et al., 2016). We calculated FLII and HFP values in QGIS for every camera trap location. We also included a random-intercept effect in our abundance formula to account for unmodeled variation between landscapes and to account for three landscapes with repeated sampling. Finally, we included a fixed effect for sampling effort (in trap nights) per sampling unit in our detection probability formula to account for multiple cameras being resampled into the same sampling unit and unequal deployment lengths. We standardized all numeric covariates (mean = 0, SD = 1) and examined Pearson correlation coefficients, and ensured no collinearity among covariates (|r| < 0.5). The abundance formula modeled the expected count relative to environmental covariates for species i at sampling unit j as follows: The abundance model for the subordinate included an additional parameter estimating the effect of the latent abundance of the dominant species: δ * N dom,j . Our detection probability formula was the same for both predator and prey and modeled the detection probability of species i at sampling unit j on sampling occasion k as: Note: We started with the original two-species N-mixture model developed by Brodie et al. (2018) (Poisson). We compared this to models using either the informed zero-inflated Poisson distribution (iZIP), or a random effect per sampling unit and sampling occasion (OD), and to the models presented in the main text using both (iZIP + OD). The species interaction column shows the posterior distribution's mean and standard deviation of the species interaction parameter, and Rhat values between 1 and 1.2 denote that the parameter successfully converged. Bayesian p-values between .25 and .75 indicate a suitable goodness-of-fit, a value of .5 indicates a perfect fit, and values outside this range indicate unacceptable performance. C-hat values greater than 1.1 indicate unacceptable overdispersion.

| Model performance
Our co-abundance modeling approach that accounted for zero inflation and overdispersion in count histories using both the iZIP distribution and ODRE successfully converged across all species pairs (Rhat <1.2; Table 2; Figure 2). Including either the iZIP and ODRE parameters improved model goodness-of-fit (Bayesian pvalues were closer to .5 than models lacking these parameters,  (Table 2b).
Co-abundance models using a Poisson instead of a binomial detection process produced lower absolute population size estimates for most species but followed the same relative abundance trends between landscapes ( Figure S1) and equivalent directionality for our species interaction parameter (δ * N dom,j ; Figure S2). However, using a Poisson detection formula models exhibited unaccep-  Figure S3) and reduced precision of the species interaction parameter that undermined its biological and applicable utility.
Therefore, we present co-abundance models using the binomial detection formula in the main text and the Poisson detection formula in the supplement.

| Southeast Asian predator-prey dynamics
We collected 5980 independent captures of our four study species over a sampling effort of 58,071 trap nights (Tables S1 and S2). We observed both clear (i.e., the 95% CI does not include 0) negative and positive predator-prey relationships from our F I G U R E 2 Comparing the effect sizes, 95% Bayesian credibility interval, and parameter convergence (Rhat values) of the species interaction parameter in all four species pairs from our two-species N-mixture models. We examined the original two-species N-mixture model proposed by Brodie et al. (2018) (Poisson), then examined the addition of the iZIP parameter (ZIP) and ODRE parameter (Poisson + OD) separately, and finally compared our final models that include both the iZIP and ODRE parameters (ZIP + OD). The Y-axis represents the mean effect sizes of the species interaction parameter, and the error bars represent the 95% Bayesian credibility interval. Finally, the values at the bottom of the bar graphs represent the specific Rhat scores from the species interaction parameter associated with each model. co-abundance models (Figures 4 and 5). The detection probability of all species assessed across all species pairs showed a clear positive relationship with sampling effort (>99% probability; Figure S4).

| Muntjac ~ tiger
The abundance of tigers clearly and negatively influenced muntjac

| Sambar ~ tiger
There was no clear relationship between tigers and sambar deer

F I G U R E 3
Comparing the goodness-of-fit between models by inspecting Bayesian p-values (a) and the magnitude of overdispersion C-hat values (b) across four species pairs. We examined the original two-species N-mixture model proposed by Brodie et al. (2018) (Poisson), then examined the addition of the iZIP parameter (ZIP) and ODRE parameter (Poisson + OD) separately, and finally compared our final models that include both the iZIP and ODRE parameters (ZIP + OD). (a) Bayesian p-values are calculated by taking the mean value of the number of times data simulated from the joint posterior distribution was greater than the real data supplied to the model, where Bayesian p-values between .25 and .75 indicate good fit, a value of 0.5 indicates a perfect fit, and values above or below the dashed red lines (<0.25 or >0.75) indicate a lack of fit. (b) C-hat values are calculated by dividing the observed data supplied to the model from data simulated from the joint posterior distribution and we visualized the mean value, where C-hat values greater than 1.1 indicate remaining overdispersion and values close to 1 indicate no remaining overdispersion. A horizontal line is added at 1 to indicate the ideal value for our C-hat scores, while the red dashed line at 1.1 denotes our cut-off point for C-hat values that suggest overdispersion. Both dominant (darker tan color) and subordinate (lighter tan color) species are included in both figures.

F I G U R E 4
Parameters describing prey species abundance from the Bayesian co-abundance models using the informed zero-inflated Poisson (iZIP) distribution and ODRE. Plots show the posterior mean effect size, and the error bars represent the 95% Bayesian credibility interval (CI), with asterisks (*) denoting relationships where the 95% CI does not include zero. The variables are FLII (green), HFP (yellow), and the species interaction (red), which shows the effect of the dominant (predator) on the subordinate (prey).

| DISCUSS ION
We introduced a two-species N-mixture modeling approach to quantify predator-prey relationships from observational count histories, while accounting for shared responses to the environment and propagating uncertainty throughout the modeling process. We assessed how this model performed using camera trapping data of two Asian apex predators and two of their key prey species that exhibit contrasting detection probabilities and natural densities. We found that explicitly classifying both the source of true zeros with the iZIP distribution in the abundance formula and false zeros with our detection formula containing the ODRE were necessary to infer ecologically meaningful predator-prey interactions while ensuring parameter convergence across all species pairs (Blasco-Moreno et al., 2019;Martin et al., 2005). Failing to classify true zeros, such as at landscapes where a species has been extirpated, leaves the model to estimate nonzero abundance due to imperfect detection (i.e., false zeros) that leads to overconfidence in the posterior estimates and increases the risk of a type I error (Martin et al., 2005). We illustrated this type I error by showing a spurious clear positive relationship between tigers and sambar deer when excluding the iZIP parameter, which disappeared when using the iZIP distribution. Our approach also produced relationships with environmental covariates across all species that are supported by past research suggesting anthropogenic pressure suppresses ungulate abundances and that intact tropical forests support large carnivore abundances Luskin, Albert, & Tobler, 2017;Macdonald et al., 2019;Rayan & Linkie, 2020).
There is a strong appetite to infer species interactions from camera trap data and to overcome the limitations of traditional cooccurrence models (Blanchet et al., 2020). While our approach is fit-for-purpose, there are several drawbacks including these models being data hungry, computationally demanding, and statistically complicated. To overcome data limitations, such as for tigers in our datasets that were only detected only in three of ten sampled landscapes, collaborative projects may be required to achieve sufficient sample sizes and span a wide gradient of predator abundances. Another solution may be to integrate multiple data sources such as direct observations with camera trapping (Miller et al., 2019). The computational power required to run and test several models with 1,000,000 MCMC iterations may pose a considerable barrier for many field ecologists (Pettorelli et al., 2021).
Access to high-performance research computing that simultaneously solves multiple models while distributing MCMC chains across multiple cores will save researchers a substantial amount of time in such analyses (Visser et al., 2015). Finally, there are concerns about the reliability of N-mixture models when assumptions are violated, such as when the same individual is observed multiple times in a single sampling occasion (Link et al., 2018).
This violation can cause inflation in absolute population size estimates, so if researchers require accurate latent population size estimates, we encourage the use of a Poisson detection process (Nakashima, 2020). We implemented co-abundance models with both binomial and Poisson detection processes and observed equivalent directionality and similar species interaction parameters, suggesting either detection process may be suitable for examining the spatial variation in abundance (Gilbert et al., 2021).
Predator losses in temperate ecosystems often lead to trophic release, where a subset of prey species increase in abundance (Ripple et al., 2009), but predator-prey relationships appear weak or positive in tropical forests (Brodie & Giordano, 2013). Predator-prey relationships in tropical forests may be strongly shaped by diffuse food webs with functionally redundant links, strong bottom-up control, and/or the joint suppression of many animals by overwhelming disturbances such as defaunation (Benítez-López et al., 2019;Brodie & Giordano, 2013;Polis & Strong, 1996;Wright et al., 1994). In accordance with this work, we found only one clear negative relationship suggesting tigers regulate muntjac deer. Further, we detected positive predator-prey relationships between clouded leopards and both muntjac and sambar deer, suggesting predation is not the dominant force regulating Asia's deer. Instead, prey abundance may be a fundamental determinant of predator abundance (Carbone & Gittleman, 2002;Karanth et al., 2004) with deer being noted as key prey for clouded leopards (Can et al., 2020;Petersen et al., 2020).
These results also suggest that hunting deer may suppress clouded leopards via prey depletion (Wolf & Ripple, 2016), or hunting may jointly suppress clouded leopards and deer together (Benítez-López et al., 2019;Ford & Goheen, 2015). Future research could apply our modeling approach to reverse the role of dominant and subordinate species to examine whether the abundance of clouded leopards is positively associated with deer abundance, thereby testing support for bottom-up regulation.
This modeling framework is widely applicable across species, ecosystems, and sampling approaches (acoustic monitoring, camera trapping, and point counts) that produce count histories. A key improvement of our co-abundance modeling approach that classifies the source of zeros in count histories is to move trophic cascades research past binary comparisons of landscapes based on predator presence or absence (e.g., Atkins et al., 2019;Brashares et al., 2010) to a gradient of predator abundance, including where predators are absent. Our modeling framework incorporates covariates and has the potential to explain divergent co-abundance patterns across environmental gradients, such as edge-adapted and disturbancetolerant species proliferating in predator-free forest fragments (Terborgh et al., 2001). Future research could use this co-abundance model to examine how predator-prey relationships vary across ecological gradients, which has important conservation implications for predicting trophic cascades (Terborgh, 2015).

| CON CLUS ION
Our study provides ecologists with a modeling approach to infer predatory interactions from observationally collected count data across multiple landscapes. This remarkably flexible approach opens novel opportunities to evaluate species interactions across natural or disturbance gradients, including when species abundances vary or where some sites experienced extirpations. There is direct applicability to ongoing predator conservation and trophic cascades research as our study species continue to experience range contractions. We look forward to opportunities to ground truth the results from our approach by co-locating camera trapping studies with manipulative experiments (e.g., predator reintroductions) and direct predation observations.

ACK N OWLED G M ENTS
In kind support was provided by the Fauna and Flora International- Animals and is copyrighted. Conversations with Jedediah Brodie, Ken Kellner, and J. Andrew Royle provided helpful insights at the early stages of this project. We thank the members of the Ecological Cascades Lab, Andrew Letten, Diana Fisher, and 2 anonymous reviewers for comments that improved previous drafts.

FU N D I N G I N FO R M ATI O N
The research was funded by the Smithsonian Institution's ForestGEO program, the Nanyang Technological University in Singapore, the University of Queensland, the National Geographic Society Committee for the Research and Exploration #9384-13, and numerous small grants for field work. M.S.L. was supported by an Australian Research Council Discovery Early Career Award no. DE210101440.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest pertaining to the research conducted here.

O PE N R E S E A RCH BA D G E S
This article has earned Open Data and Open Materials badges.

DATA AVA I L A B I L I T Y S TAT E M E N T
All R code and data used to implement our Bayesian co-abundance models can be accessed at this Figshare repository: https://doi. org/10.5061/dryad.b8gtht7h3.